Phase diagram of the two-chain Hubbard model 
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We have calculated the charge gap and spin gap for the two-chain Hubbard model as a function 
of the on-site Coulomb interaction and the interchain hopping amplitude. We used the density 
matrix renormalization group method and developed a method to calculate separately the gaps 
numerically for the symmetric and antisymmetric modes with respect to the exchange of the chain 
indices. We have found very different behaviors for the weak and strong interaction cases. Our 
calculated phase diagram is compared to the one obtained by Balents and Fisher using the weak 
coupling renormalization group technique. 
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Although the Luttinger liquid behavior of the one- 
dimensional Hubbard model has been well understood, 
the two-dimensional Hubbard model, which is believed 
to be related to understanding of the high T c supercon- 
ductivity jjj, is not yet clear. As a crossover between ID 
and 2D, the two-chain Hubbard model is certainly a good 
theoretical basis for the ladder compounds . It is im- 
portant to understand how two Luttinger liquid systems 
evolve to the ladder system as the interchain coupling is 
introduced. 

The Hamiltonian for the two-chain Hubbard model is 



H = ~t\\ ^ ( C L°y<7 + H - C ) - f -L ^2( C \i<r C 2ia + H.C 

l,{i,j),a i,(r 

+ U ^2 n U,l n U,U 



(1) 



where I is the chain index, I = 1,2. We have the intra- 
chain hopping t\\ term, the interchain hopping t± term 
and the on-site Coulomb interaction U term, tu = 1 in 
this paper. 

Some authors have studied the phase diagram of the 
two-chain Hubbard model using the weak coupling renor- 
malization group method f§-^|. At half-filling, the sys- 
tem is an insulator with a spin gap. Upon light hole 
doping, for small t± the spin gap remains finite. For 
large t±, the complete separation between the bonding 
band and the antibonding band leads the system to a 
Luttinger liquid phase. But it remains unclear for large 
U case which is equivalent to the t — J ladder @-|l(| . Ba- 
lents and Fisher described the phase in terms of number 
of gapless charge and spin modes, which they denoted 
by CnSm where n is the number of gapless charge modes 
and m is the number of gapless spin modes Q . For our 
work we also use this notation. Their phase diagram is 
rather diverse in the hole-doped region. Interestingly, 
they found phases such as C2S2 and C2S1 between the 
C1S1 Luttinger liquid phase and the spin-gapped CI SO 



phase. And the phase diagram is qualitatively same for 
different electron filling n except for the cases of half- 
filling, quater- filling and half-filled bonding band, where 
the umklapp process will be relevant. 

Since most of the previous results done by the weak 
coupling renormalization group method may be reliable 
only in the U — > + limit, the questions are whether it 
will be similar for finite U and whether it will depend on 
U. Noack et al. have studied the correlation functions 
and the gaps for finite U using the density matrix renor- 
malization group method (DMRG) ]ll| . They found the 
enhancement of the d-wave pairing correlation function 
for the spin-gapped phase at rather anisotropic regime 
t± > in, though they cannot distinguish the symmetric 
and antisymmetric modes. Hence it is difficult to deter- 
mine which CnSm phase it is. To understand the phase 
diagram in detail and to make comparison with the weak 
coupling result it is necessary to study the charge and 
spin excitations in different modes for finite U. 

In this paper, we study the phase diagram of the two- 
chain Hubbard model by calculations of the charge and 
spin gaps for the symmetric mode and antisymmetric 
mode respectively for finite values of U at a fixed electron 
filling ri = 0.75. Using this result we can easily under- 
stand the result for other electron filling. For our study, 
we use the density matrix renormalization group method 
(DMRG) JiJ,[l|] and develop a method to differentiate 
the modes. 

We take the transform of c operators in terms of bond- 
ing and antibonding forms as following: 
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Then the Hamiltonian becomes 
H =-t\\ ( c Ux c +^ 

(id)," 



1 



-*L ^2(cl ia c +la - cl la c- la ) 

+j Yl^ n +^ + n -n)( n +n + n -a) 
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+ c + J T c + J i c -a c -n + c^ iT cL a c +ii c +iT 
+ c W-a c +a c -*T + cL^c^c-^c+it]. (3) 

Besides the usual symmetry of the total spin S and S z , 
the translational symmetry along the chain direction, the 
Hamiltonian also has the symmetry of the exchange of 
the chain indices. The symmetry of the total state under 
the exchange of two chain indices depends on the sym- 
metries of wavefunction of each particle. It is decided 
by the number of particles on the odd chain of the new 
transformed Hamiltonian, since the evenness and odd- 
ness of the number of particles on the odd chain is not 
changed by the Hamiltonian. If the number of particles 
on the odd chain is odd, total state is antisymmetric and 
if the number is even, the total state is symmetric. In 
the DMRG calculation, we define the symmetry under 
the chain exchange of each base or sector in terms of 
the number of particles on the odd chain and so we have 
the option to choose the symmetry and able to calculate 
the ground state energy of a given symmetry. We also 
use the open boundary condition here hence there is no 
translational symmetry along the chain direction. 

For a set of given values of U, tj_ and electron filling 
n, the charge gaps and spin gaps are defined as follows: 

A C+ = ±[E+(Q + 2,S = S Z = 0) 
+ E+(Q-2,S = S Z =0)- 2E+(Q, S = S z =0)], (4) 



mode. But for S z = sector with antisymmetric mode, 
the ground state is S = 1 state rather than S = 0, which 
implies A c _ > A s _ for finite size lattices. Therefore, to 
calculate the antisymmetric charge gap we need to find 
the S = state which is an excited state. We find this 
state by calculating the expectation values of S for a few 
excited states. 
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FIG. 1. Gaps versus 1/L to calculate the thermodynamic 
values of gaps for U = 8 and t± = 1.7. We denote 
symmetric charge gaps with dark circles, antisymmetric 
charge gaps with diamonds, symmetric spin gaps with 
open circles and antisymmetric spin gaps with squares. 



A c _ = E_(Q, S = S z = 0) -E + (Q, S = S Z = 0), (5) 
A s+ - E+(Q, S = S Z = 1)- E+{Q, S = S Z = 0), (6) 

A s _ = E_(Q, S = S z = l)-E+(Q,S = S z = 0), (7) 

where ± is symmetric or antisymmetric mode and all en- 
ergies are the lowest energy for each set of quantum num- 
bers. Here, we assume the total number of the particles is 
even. In each iteration of the DMRG calculation, we typ- 
ically keep M = 200 states for the block. We calculated 
gaps for various size of lattices with length L = 8, L = 12, 
L = 16 and L = 20 and extrapolated the thermodynamic 
values by 1/L polynomial expansions. In the DMRG cal- 
culation, we only utilize the conservation of S z , not the 
total spin S. To identify the total spin S correctly, it is 
necessary to check S values of the lowest energy state in 
each sector. For S z = sector with symmetric mode, 
the ground state is S = state. For S z = 1 sector with 
symmetric mode, the ground state is S = 1 state and 
the same is true for S z = 1 sector with antisymmetric 



Fig. 1 shows the finite size extrapolation used to ob- 
tain the thermodynamic values of the gaps for U = 8 
and tj_ = 1.7. The symmetric charge gap A c+ and the 
symmetric spin gap A s+ vanish. Specially, for the entire 
range of U and t±, A c+ always vanishes at n = 0.75. 
The antisymmetric charge gap A c _ and the antisymmet- 
ric spin gap A s _ are finite. A c _ is larger than A s _ for 
all size of lattices but, for the longer lattice the difference 
between them is smaller. Since we take the ground state 
as the target state for the calculations of the density ma- 
trix in DMRG, the accuracy of the excited states is not 
as high as that of the ground state and this is even more 
serious for longer lattice and for small t±. 

For strong interaction U = 8, A s+ and A s _ are plotted 
in Fig. 2. as functions of t± for various size of lattices 
including L — oo extrapolated values. For most values 
of t± we get very good extrapolations for both A s+ and 
A s _with three L values: 8, 12 and 16. There are some 
points near t± = 1.0 for which we are not able to get a 
good extrapolation. As t± decreases from t± = 2.0 A s _ 
decreases but does not vanish until t± = 0.6, and A s+ 
opens up around at t± = 1.5. Since we always have one 
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gapless symmetric charge mode, that is A c+ = and 
A c _ > A s _ the phase changes from C1S1 phase to C1S0 
phase around t± = 1.5. Between these two phase we do 
not have phases such as C2S2 and C2S1 which are found 
in the results for U — > + limit by Balents and Fisher || . 
We stop our calculations at t± = 0.6 because the DMRG 
calculation has poor accuracy for small t±. 
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which is the Luttinger phase, and we will have C2S2 
phase with all gapless modes. So, as t± decreases all gaps 
will vanish one by one and eventually we will have C2S2 
phase in the t± — limit. The linear behavior of A s _ in 
the C1S1 phase region for t± > 1.5 can be understood in 
the limit of infinitely large t±. In this limit, A s _ mea- 
sures the separation between bonding band and empty 
antibonding band. In the physically interesting isotropic 
region t± = 1.0, both spin gaps are similar in magni- 
tude but finite, which is consistent with previous results 
pj . Also, they increase slowly as t± increases. This is 
compatible with the results for the t — J ladder, which is 
linear behavior of the spin gap in J± ( A s ~ J± ~ t\ /U 
) [f|J§l. The increasing behavior of A s+ for t± < 0.7 is an 
artifact from poor convergence of DMRG for small t± . 
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FIG. 2. Spin gaps versus t± for U = 8. (a) for symmetric 
mode and (b) for antisymmetric mode. We denote values 
for L = 8 with diamonds, L = 12 with squares, L = 16 
with open circles and L — oo extrapolated values with 
dark circles. 

The vanishing of A s _ at t± = 0.6 can be understood 
from the t± = limit. In this limit, two chains will 
be completely decoupled to two one-dimensional chains, 
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FIG. 3. Same figures as FIG.2 for U = 1. 
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Fig. 3 shows same figures for weak interaction U = 1 
as Fig. 2. Unlike U — 8 case, A s+ does not open up 
and A s _ vanishes around at t± = 1.7, where A c _ also 
vanishes so that we have a transition from C1S1 phase 
to C2S2 phase. As t± decreases further, we have finite 
size fluctuations so that we can not calculate L = oo 
gaps until tj_ = 0.8. Surprisingly, in a small region of 
t±, 0.6 < t± < 0.8, there is little finite size fluctuation. 
We find both A s+ and A s _ vanish. On the other hand 
A c _ is finite, but considering the overestimation of A c _ 
because of the inaccuracy of the excited states and the 
system being two decoupled Luttinger liquid systems for 
small t±, presumably it also vanishes and we have C2S2 
phase. If we smoothly connect this region to the region 
with t± ~ 1.7, it is plausible to assume that the phase 
remains gapless. Therefore, in this weak interaction case, 
there will be only one transition from C1S1 phase at large 
t± to C2S2 phase at small t± and this is rather consistent 
with the nonintcracting limit. In the C1S1 phase, we also 
have the linear behavior of A s _ with smaller A s _ than 
U = 8 case but with larger slope. 



When we have different n, the transition point of t± will 
be changed because it is around the value of the nonin- 
teracting limit, t± = 1 — cos7rn just as the weak coupling 
case ||. Therefore, for different n the whole phase dia- 
gram will be shifted in t± direction. 

In conclusion, we have derived the phase diagram of 
the two-chain Hubbard model in terms of charge and spin 
excitations for both symmetric and antisymmetric modes 
at fixed electron filling. We have found very different 
behaviors for the weak and strong interaction cases. The 
transition from the noninteracting limit to the interaction 
case seems to be rather gradual and beyond finite value of 
U we have the spin-gapped phase at the isotropic region 
t± = 1.0, which is consistent with the previous results 
for the strong interaction case and t — J model. 
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Rep. of China, Grant Nos. NSC87-2112-M-011-016. 
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FIG. 4. Phase diagram at n = 0.75. We denote C1S1 
phase with open circles, C1S0 with squares, C2S2 with 
dark circles and C2S1 with diamonds. We also show the 
phases for the noninteracting case. 

With calculations for other values of U, we plot the 
phase diagram in the space of t± and U at n = 0.75 
in Fig. 4. It is clear that the strong interaction case 
(U = 8) and the weak interaction case (U — 1) are 
clearly different. {7 = 2 case is same as U — 1 case 
and U = 4 seems to be similar to U = 8 case regardless 
of the existence of the finite size fluctuations. The tran- 
sition point of i_L between two different phases decreases 
from t± ~ 1.7 for the noninteracting case as U increases. 



46, 



50, 



12133 
Phys. 



[1] P.W. Anderson, "The theory of Superconductivity in the 
High-T c Cuprates" (Princeton University Press, Prince- 
ton, 1997). 

[2] E. Dagotto and T.M. Rice, Science 271, 618 (1996). 
[3] D.J. Scalapino, Physica C 282, 157 (1997). 
[4] M. Fabrizio, A. Parola and E. Tosatti, Phys. Rev. B. 
3159 (1992). 

[5] D.V. Khveshchenko and T.M. Rice, Phys. Rev. B. 
252 (1994). 

[6] L. Balents and M.P.A. Fisher, Phys. Rev. B 53, 
(1996); H.H. Lin, L. Balents and M.P.A. Fisher, 
Rev. B 56, 6569 (1997). 

[7] E. Dagotto, J. Riera and D. Scalapino, Phys. Rev. B 45, 
5744 (1992). 

[8] M. Sigrist, T.M. Rice and F.C. Zhang, Phys. Rev. B 49, 
12058 (1994). 

[9] H. Tsunetsugu, M. Troyer and T.M. Rice, Phys. Rev. B 

49, 16078 (1994). 
[10] D. Poilblanc, H. Tsunetsugu and T.M. Rice, Phys. Rev. 

B 50, 6511 (1994); D. Poilblanc, D.J. Scalapino and W. 

Hanke, Phys. Rev. B 52, 6796 (1995). 
[11] R.M. Noack, S.R. White and D.J. Scalapino, Physica C 

270, 281 (1996); R.M. Noack, N. Bulut, D.J. Scalapino 

and M.G. Zacher, Phys. Rev. B 56, 7162 (1997). 
[12] S.R. White, Phys. Rev. Lett. 69, 2863 (1992); Phys. Rev. 

B 48, 10345 (1993). 
[13] S. Liang and H. Pang, Phys. Rev. B 49, 9214 (1994); 51, 

10287 (1995). 



4 



